#ifndef __MISC_PROCESSING_HPP
#define __MISC_PROCESSING_HPP

// misc_processing.hpp
// 包含各种用于数据处理的函数

#include <opencv2/core.hpp>
#include <algorithm>
#include <valarray>
#include <iterator>

// 一维离散傅里叶变换
template <typename IT>
cv::Mat opencv_fft1(IT begin, IT end);

// 在指定范围寻找峰值
template <typename IT>
int find_peak(IT begin, IT end);

// 计算复数矩阵各个元素的模
cv::Mat complex_abs(cv::Mat complex_mat);

// 计算标准差
template <typename IT>
double std_dev(IT begin, IT end);

// 将数据正则化为最大值为y_max，最小值为y_min
template <typename IT, typename T>
std::valarray<T> normalized(IT begin, IT end, T y_min, T y_max);

// 计算图像某一区域颜色的均值
cv::Vec3d area_avg(cv::Mat mat, cv::Rect area);

// 色差色度法（build chrominance signals）
// 输入应该是指向cv::Vec3x的迭代器
template <typename IT>
std::valarray<double> chrominance(IT begin, IT end, double k_r, double k_g, double k_b);

// 计算平均值
template <typename T, typename IT>
T average(IT begin, IT end);

template <typename IT>
cv::Mat opencv_fft1(IT begin, IT end) {
    int size = end - begin;
    cv::Mat ans(1, size, CV_64F);
    IT p;
    cv::MatIterator_<double> q;
    for (p=begin,q=ans.begin<double>(); p!=end; ++p,++q) {
        *q = *p;
    }

    cv::Mat planes[] = {cv::Mat_<double>(ans), cv::Mat::zeros(ans.size(), CV_64F)};
    merge(planes, 2, ans); 

    cv::dft(ans, ans);

    return ans;
}

template <typename IT>
int find_peak(IT begin, IT end) {
    int ans = 0;
    IT p, q;
    for (p=begin, q=begin+1; q!=end && *p<*q; ++p, ++q) {
        ++ans;
    }
    return ans;
}

template <typename IT>
double std_dev(IT begin, IT end) {
    double sqr_avg = 0.0; // 平方的均值
    double avg_sqr = 0.0; // 均值的平方
    int n = end -begin;

    for (IT p=begin; p!=end; ++p) {
        avg_sqr += *p;
        sqr_avg += (*p)*(*p);
    }
    avg_sqr /= n;
    avg_sqr *= avg_sqr;
    sqr_avg /= n;

    return std::sqrt(sqr_avg-avg_sqr);
}

template <typename IT, typename T>
std::valarray<T> normalized(IT begin, IT end, T y_min, T y_max) {
    T x_max = *std::max_element(begin, end);
    T x_min = *std::min_element(begin, end);
    T delta_x = x_max - x_min;
    T delta_y = y_max - y_min;
    int n = end - begin;

    std::valarray<T> ans;
    ans.resize(n);

    IT p;
    int i;
    for (i=0, p=begin; i<n; ++i, ++p) {
        ans[i] = delta_y/delta_x * (*p-x_min) + y_min;
    }

    return ans;
}

template <typename IT>
std::valarray<double> chrominance(IT begin, IT end, double k_r, double k_g, double k_b) {
    int n = end - begin;
    std::valarray<double> ans;
    ans.resize(n);

    int i;
    IT p;
    for (i=0, p=begin; i<n; ++i, ++p) {
        ans[i] = k_b*(*p)[0] + k_g*(*p)[1] + k_r*(*p)[2];
    }

    return ans;
}


template <typename T, typename IT>
T average(IT begin, IT end) {
    T sum = 0;
    for (IT p=begin; p<end; ++p) {
        sum += *p;
    }
    return sum / (end - begin);
}

#endif // __MISC_PROCESSING_HPP